Entrainment of randomly coupled oscillator networks by a pacemaker 
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Entrainment by a pacemaker, representing an element with a higher frequency, is numerically 
investigated for several classes of random networks which consist of identical phase oscillators. We 
find that the entrainment frequency window of a network decreases exponentially with its depth, 
defined as the mean forward distance of the elements from the pacemaker. Effectively, only shallow 
networks can thus exhibit frequency-locking to the pacemaker. The exponential dependence is also 
derived analytically as an approximation for large random asymmetric networks. 
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Pacemakers are wave sources in distributed oscillatory 
systems typically associated with a local group of ele- 
ments having a higher oscillation frequency. Target pat- 
terns, generated by pacemakers, were the first complex 
wave patterns observed in the Belousov-Zhabotinsky sys- 
tem |]|. Pacemakers play an important role in func- 
tioning of the heart |2j and in the collective behavior 
of Dictyostelium discoideum Q. They are also observed 
in large-scale ecosystems Q. In addition to pacemak- 
ers produced by local heterogeneities in the medium ||, 
self-organized pacemakers in uniform birhythmic media 
have been theoretically studied [||. While the majority of 
related investigations have so far been performed for sys- 
tems with local diffusive coupling between the elements, 
pacemakers can also operate in oscillator networks with 
complex connection topologies. The circadian rhythm 
in mammals is a daily variation of 24 hours that reg- 
ulates basic physiological processes in such animals []|. 
It is produced by a complex network of neurons form- 
ing the so-called suprachiasmatic nucleus (SCN) As 
recently shown, this oscillator network undergoes spon- 
taneous synchronization in absence of any environmen- 
tal input, but its intrinsic synchronization period is then 
significantly longer than 24 hours 0. Therefore, the ac- 
tual shorter rhythm results from the environmental en- 
trainment and must be externally imposed. The entrain- 
ment is mediated by direct photic inputs from eyes into 
the SCN, which undergo periodic daily variation. How- 
ever, it is known that only a distinct subset of neurons 
in this network is directly influenced by photic inputs 
|lfj| . Hence, functioning of this particular neural system 
is crucially dependent on the ability of the entire complex 
network to become entrained by an external pacemaker. 
Analogous behavior can also be expected, for example, in 
heterogeneous arrays of globally coupled electrochemical 
oscillators where synchronization and entrainment have 
been experimentally demonstrated |TT| . 

To understand operation of pacemakers in networks 
with complex connection topologies, action of a pace- 
maker in a random oscillator network should first be in- 
vestigated. In this Letter, networks of identical phase 



oscillators with random connections are considered. A 
pacemaker is introduced as a special element whose os- 
cillations have a higher frequency and are not influenced 
by the rest of the system. Depending on the pacemaker 
frequency and the strength of coupling, the pacemaker 
can entrain the entire network, so that the frequencies 
of all its elements become equal to that of the pace- 
maker. We find that the entrainment window decreases 
exponentially with the depth of a network, defined as 
the mean forward distance of its elements from a pace- 
maker, and thus only shallow networks can effectively be 
entrained. This result is confirmed in numerical simu- 
lations for several different classes of random networks, 
including small-world graphs. It is further analytically 
derived as an approximation for random networks with 
asymmetric connections. 

We consider a system of N + 1 phase oscillators, one of 
them being a pacemaker. The model is given by a set of 
evolution equations for the oscillator phases fa and 
the pacemaker phase fa, 
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sin(fa - fa) - /sBi sm(fa -fa), 

(1) 



The topology of network connections is determined by 
the adjacency matrix A whose elements Aij are either 1 
or 0. The element with i = is special and represents 
a pacemaker. Its frequency is increased by Auj with re- 
spect to the frequency u> of all other oscillators 0] . The 
pacemaker is acting on a randomly chosen subset of iVi 
elements, specified by Bi taking values 1 or 0. The total 
number of connections to the pacemaker, N± — $^-B», 
is fixed. The coupling between elements inside the net- 
work is characterized by strength k. The strength of 
coupling from the pacemaker to the network elements is 
determined by the parameter /i. In absence of a pace- 
maker, such networks undergo autonomous phase syn- 
chronization at the natural frequency lu. Without loss 
of generality, we put ui — 0. Moreover, we rescale time 
as t' = t Alu and introduce rescaled coupling strengths 
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FIG. 1: (color online) Phases of elements in the entrained 
state; N = 100, p = 0.05, Ni = 3, k = 100 and /i = 1000. 

k? = n/Aui,fi' = fi/Au). After such rescaling, the model 
takes the form of Eqs. (1) with Auj = 1 and uj = (we 
drop primes in the notations for the rescaled couplings). 
In terms of the original model (1), increasing the rescaled 
coupling between the elements is equivalent either to an 
increase of coupling n or to a decrease of the relative 
pacemaker frequency Aco. 

The presence of a pacemaker imposes hierarchical or- 
ganization. For any node i, its distance h with respect 
to the pacemaker is given by the length of the minimum 
forward path separating this node from the pacemaker. 
All Ni elements in the group directly connected to the 
pacemaker have distances h = 1, the next elements which 
are connected to the elements from this group have dis- 
tances h — 2, etc. Thus, the whole network is divided 
into a set of shells [l^, each characterized by a certain 
forward distance h from the pacemaker. The set of num- 
bers Nh is an important property of a network. The 
depth L of a given network, which is the mean distance 
from the pacemaker to the entire network, is introduced 
as L = (l/N)J2h.hNh- It should be noticed that such 
ordering of network nodes is based solely on the forward 
connections down the hierarchy and does not depend on 
the distribution of reverse (upward) connections in the 
system. 

First, we investigated standard random asymmetric 
networks, where independently for all connections A^ = 
1 with probability p and Ajj = otherwise. Only sparse 
random networks with relatively low mean connectivity p 
and a small number Ni of elements directly connected to 
the pacemaker were considered. Numerical simulations 
were performed for the networks of size N = 100 starting 
with random initial conditions for the phases of all oscilla- 
tors. For each oscillator, its effective long-time frequency 
u)i was computed as u>i — T^ 1 [<fti(to + T) — c/>i(io)] with 
sufficiently large T and to. The simulations show that 
the response of a network to the introduction of a pace- 
maker depends on the strength k of coupling between the 



oscillators. When this coupling is sufficiently large (and 
coupling fi to the pacemaker is also sufficiently strong as 
assumed below), the pacemaker entrains the whole net- 
work (i.e., u>i — 1 for all elements i). The frozen relative 
phases ipi = <fii — 4>o are displayed in Fig. 1. Here, the 
elements are sorted according to their hierarchical shells. 
Despite random variations, there is a clear correlation be- 
tween phases of oscillators and their positions in the hi- 
erarchy. Generally, the phase decreases for deeper shells, 
and the phase difference between the neighboring shells 
rapidly becomes smaller as deeper shells are considered. 
As the coupling strength n is decreased, the entrainment 
breaks down at a certain threshold value n cr . Our simula- 
tions show that synchronization between the first and the 
second shells was almost always the first to break down, 
and the frequencies of the second and deeper shells re- 
mained equal in most cases for the considered random 
networks. 

Figure 2 displays in the logarithmic scale the thresh- 
olds k ci for a large set of networks with different depths 
and different numbers of elements in the first shell. Each 
group with a certain Ni is displayed by using its own 
symbol. Every such group generates a cluster of data 
points. Correlation between the entrainment threshold 
and the network depth is apparent. The distributions 
inside each cluster and the accumulation of the clusters 
yield the dependence K cr (L) of the entrainment threshold 
on the network depth. Note that the statistical varia- 
tion of the data becomes larger for deeper networks with 
larger L and for smaller pN. Similar dependence was 
found for the networks with different mean connectiv- 
ity p (see inset). Remarkably, the observed dependences 
could be well numerically approximated by the exponen- 
tial dependence 

K CI cx (1+ P N) L . (2) 

As the second class, asymmetric small-world networks 
fl4| were considered. To generate them, we first con- 
structed a one-dimensional lattice of N elements where 
each element had incoming connections from up to its 
fcth neighbor (the degree was thus 2k). Then a randomly 
chosen link in the lattice was eliminated and a distant 
connection between two independently randomly chosen 
elements was introduced. This construction was repeated 
qN times, with the parameter q specifying the random- 
ness of a network. When q was small, the network was 
close to a lattice and, in this case, we have seen that sta- 
ble wave solutions with different winding numbers were 
possible, depending on initial conditions (cf. 0,0). To 
avoid this, we chose almost synchronized states as initial 
conditions. The entrainment thresholds for such small- 
world networks are displayed in Fig. 3 and again show a 
clear correlation between k ci and L. The dependence on 
the depth is approximately linear in lattices (q = 0), but 
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FIG. 2: (color online) Dependence of the entrainment thresh- 
old on the depth L for an ensemble of random networks with 
N = 100 and p = 0.1. In the inset, respective data for net- 
works with p = 0.06 and 0.2 is plotted. Solid lines are the 
exponential functions c p (l + pN) L with appropriate fitting 
parameters c p . 
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FIG. 3: (color online) Dependence of the entrainment thresh- 
old on the depth L for an ensemble of small-world networks 
with N = 100, k — 3 and different randomness q. The solid 
line is the same exponential function as that fitted to the data 
with pN — 6(= 2k) in Fig. 2. The dotted line is linear fitting 
for the regular lattice (5 = 0). 



it becomes strongly nonlinear even when small random- 
ness is introduced. For q — 0.1, the dependence is already 
approximately exponential, though the dispersion of data 
is strong. As randomness q is increased, the dependence 
approaches that of the standard random networks with 
pN = 2k. 

We have also investigated asymmetric scale-free ran- 
dom networks |16(, asymmetric regular random networks 
(where every element has exactly the same number of 
either incoming or outgoing connections), and symmet- 
ric standard random networks. For all of them, ap- 
proximately exponential dependences of the entrainment 
threshold on the network depth were observed in a large 
parameter region. 

The exponential dependence @ can be approximately 
derived for asymmetric random networks with large ./V 
and pN. In the large-size limit, random graphs have 
locally a tree-like structure . The global tree approx- 
imation has previously been used for determining statis- 
tical properties of random networks |l6l | . We apply here 
the same approximation and assume that the graph of 
forward connections extending from the pacemaker node 
represents a tree, so that any oscillator has only one in- 
coming connection from the previous shell. Then the 
shell populations Nh are given by Nh = Ni(pN) h ^ 1 for 
h = 2, . . . , H, where H is the total number of shells de- 
termined by J2h=i^ h = N. Because pN is large, we 
have Nh <C Nh — N for h < H, and thus L ~ H. Next, 
we estimate the numbers mhk of incoming connections 
leading from all elements in the fcth shell to an oscilla- 
tor in the hth shell. By definition of hierarchical shells, 
fnhk — if k < h—1. In the tree approximation, mhk = 1 



for k = h — 1. Because most of the population is concen- 
trated in the last shell, reverse connections from other 
shells can be neglected. On the average, the number of 
reverse connections from the shell H to an oscillator in 
the shell h is rrihH = pNh- Moreover, the relative sta- 
tistical deviation from this average is of order {pN)~ x / 2 , 
and is thus negligible. Therefore, in this approximation 
all oscillators inside a particular shell have effectively the 
same number of connections from other shells, and a state 
with phase synchronization inside each shell is possible. 
In this state, all oscillators inside a shell have the same 
phase, i.e. fa = Oh for all oscillators i in a shell h. Un- 
der entrainment, the phases of such a state can be found 
analytically as a solution of algebraic equations 

H 

- — — y~] m hk - sm.(6 h -0 k ) =1 for h = 2, . . . , n, 
P N k^-i 

H 

- /isin(6»i -0 O ) — y~] mik sin(0i - k ) = 1, (3) 

pN 

where $o = 4>o- For large pN , we can linearize sin ( 0^ — 0k ) 
for h, k > 2 in the solution of Eqs. 10 [it can be shown 
that O2 — Oh is of order 0(l/pN)]. Furthermore, using 
that N H o± N > N h for h < H and L ~ H, we obtain 
for ft, > 2 that 

sm(0h-i-0h) = —(l+ P N) L - h . (4) 

K 

Equations (4) determine the phases of oscillators in the 
considered synchronized state. Note that the explicit 
value of the phase 0\ in this state is not needed below. 
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The entrainment breakdown can, in principle, oc- 
cur through destabilization of the synchronized state. 
Though the analytical proof of its stability is not yet 
available, our numerical simulations show that the syn- 
chronous entrained state with < 9h-i — Oh < tt/2 is 
always stable when it exists. Thus, the breakdown of 
entrainment in the considered system takes place in a 
saddle-node bifurcation, through the disappearance of so- 
lutions of Eqs. ©■ This occurs when | sin(0fc_i — 6h)\ = 1 
for certain h. For large enough /i, we always have 
< sin(#o — #i) < 1 (a sufficient condition is fx > 1 + k). 
Among the other terms, the term sin(0i — 62) is always 
the largest one. Therefore, the solution disappears and 
breakdown occurs when sin(#i — 62) — 1. Substituting 
Eq. 10} into this equation and solving it with respect to 
k, we finally derive the dependence (J2J. Thus, we see 
that the entrainment breakdown occurs through the loss 
of frequency locking between the first shell and the rest 
of the network. The analytical estimate for the critical 
coupling strength, obtained using the tree approxima- 
tion, agrees well with the numerical data, even for the 
networks which are not very large. 

So far we have used the coupling strength which was 
rescaled as n — > k/Alj. Therefore, if the non-scaled 
coupling strength is fixed, Eq. J5J determines the max- 
imum Aw c at which the entrainment is still possible, 
Aw c oc k(1 +pN)~ L . The entrainment by a pacemaker 
can take place only if its frequency lies inside the interval 
(a;, lu + Au> c ). 

Thus, the entrainment window decreases exponentially 
with the depth of a network. This is the principal re- 
sult of our study, which holds not only for standard ran- 
dom networks, where the above analytical estimate is 
available, but also for small-world graphs and other nu- 
merically investigated random topologies. In practice, it 
implies that only shallow random networks with small 
depths are susceptible to frequency entrainment. 

Our results remain valid when, instead of a pacemaker, 
external periodic forcing acts on a subset of elements. We 
have checked that the reported strong dependence on the 
network depth remains valid for systems with larger net- 
work sizes, heterogeneity in frequencies of individual os- 
cillators, and several other coupling functions. The study 
was performed for coupled phase oscillators which serve 
as an approximation for various real oscillator systems, 
including neural networks (see, e.g., 0, 0, llflj)- Its 



conclusions should be applicable for a broad class of os- 
cillator networks with random architectures. 
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